
#download necessary libraries
import numpy as np
import datetime as dt
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
from dateutil.relativedelta import relativedelta
from scipy.fftpack import fft, fftfreq
import cartopy.crs as ccrs
import cartopy.feature as cfeature

#define style
%matplotlib inline
plt.style.use('seaborn-v0_8-colorblind')


#define the spatial reference for the map
provinces_50m = cfeature.NaturalEarthFeature('cultural',
                                             'admin_1_states_provinces_lines',
                                             '50m',
                                             facecolor='none')

land_50m = cfeature.NaturalEarthFeature('physical', 'land', '50m',
                                        edgecolor='k',
                                        facecolor=cfeature.COLORS['land'])

#create a figure and cubplot
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(111, projection=ccrs.LambertConformal())

#add an image 
ax.stock_img()

#plot the sites based on the coordinates defined in the dataset
ax.plot(-83.330, 41.762, '*',  color = 'maroon', transform=ccrs.PlateCarree(), label = "WE2")
ax.plot(-83.193, 41.827, '*', color = 'orangered', transform=ccrs.PlateCarree(), label = 'WE4')
ax.plot(-83.385, 41.705, '*', color = 'orange', transform=ccrs.PlateCarree(), label = 'WE6')
ax.plot(-83.364, 41.834, '*', color = 'gold', transform=ccrs.PlateCarree(), label = 'WE8')

ax.plot(-83.424, 41.718, '*', color = 'green', transform=ccrs.PlateCarree(), label = "WE9")
ax.plot(-83.254, 41.703, '*',  color = 'blue', transform=ccrs.PlateCarree(), label = 'WE12')
ax.plot(-83.136, 41.741, '*', color = 'blueviolet', transform=ccrs.PlateCarree(), label = 'WE13')

#add all necessary features
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.LAKES)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.RIVERS)
ax.add_feature(cfeature.BORDERS)
ax.add_feature(provinces_50m)
#add legend
ax.legend()

#set extent
ax.set_extent([-82, -85, 40, 43], ccrs.Geodetic())


#convert times to a datetime object when they are read in the format month-day-year

def tconvert(x):
     '''
    returns a datetime object for a string in the format month-day-year
    '''

     return dt.datetime.strptime(str(x), '%m/%d/%y')
    
#define a function that will read in the files
def nutrient_data_file_reader(filename):
     '''
    reads in a comma-separated file using genfromtxt that has dates in the format month-day-year which 
    will be converted to datetime objects using the converter defined above. File column names will be used to 
    parse the data
    '''

     return(np.genfromtxt(filename, delimiter = ',', converters = {0: tconvert}, encoding = 'utf-8',
                          names = True, dtype = None))
        
#create an array of all the file names for looping
filenames = ['WE2.csv', 'WE4.csv', 'WE6.csv', 'WE8.csv', 'WE9.csv', 'WE12.csv', 'WE13.csv']

#create an empty dictionary of nutrient data
nutrient_data = {}

#loop through all the files and add them to the dictionary of nutrient data
#the split function allows us to break up the filename so the key is just the site name without the .csv ending
for name in filenames:
    key = name.split('.')[0]
    nutrient_data[key] = nutrient_data_file_reader(name)


# we are interested in Nitrogen and Phosphorus concentrations at these sites so first we will plot this nutrient data
# for each of these sites
fig, (ax1, ax2) = plt.subplots(2, 1, figsize = (10, 12))
plt.subplots_adjust(hspace=0.25)
# plot each site color the same as the map
count = 0
colors = ['maroon', 'orangered', 'orange', 'gold', 'green', 'blue', 'blueviolet']
for name in filenames:
    key = name.split('.')[0]
    color = colors[count]
    count = count + 1
    ax1.plot(nutrient_data[key]['\ufeffDate'], nutrient_data[key]['Nitrate_Nitrite'], label = key, color = color)
ax1.legend(loc='upper left')
ax1.set_xlabel('Year', fontsize='x-large')
ax1.set_ylabel('N/N Concentrations [$mg N/L$]', fontsize='large')
ax1.set_title('N/N at Sites in Western Lake Erie Basin 2012-2018', fontsize = 15)

count = 0
for name in filenames:
    key = name.split('.')[0]
    color = colors[count]
    count = count + 1
    ax2.plot(nutrient_data[key]['\ufeffDate'], nutrient_data[key]['Total_Dissolved_Phosphorus'], label = key, color = color)
ax2.legend(loc='upper left')
ax2.set_xlabel('Year', fontsize='x-large')
ax2.set_ylabel('TDP Concentrations [$\mu$g P/L]', fontsize='large')
ax2.set_title('TDP at Sites in Western Lake Erie Basin 2012-2018', fontsize = 15)
plt.show()


# Plot nitrogen and phosphorus at the same site on the same plot for comparison 

#create a figure with 1 subplot
fig, (ax1) = plt.subplots(1, 1, figsize = (10, 6))

#plot the nitrogen data
ax1.plot(nutrient_data['WE6']['\ufeffDate'], nutrient_data['WE6']['Nitrate_Nitrite'], label = 'Nitrate/Nitrite',
        color = "forestgreen")

#add label
ax1.set_ylabel('N/N Concentrations [$mg N/L$]', fontsize='x-large')

#create a twin axis 
ax2 = ax1.twinx()

#plot the TDP data
ax2.plot(nutrient_data['WE6']['\ufeffDate'], nutrient_data['WE6']['Total_Dissolved_Phosphorus'],
         label = 'Phosphorus', color = "turquoise")

#add plot elements (labels, legend, title) 
ax1.set_xlabel('Year', fontsize='x-large')
ax2.set_ylabel('TDP Concentrations [$\u00B5g P/L$]', fontsize='x-large')
ax1.set_title('Nutrient Concentrations at Site WE6 in Western Lake Erie Basin 2012 - 2018', fontsize = 16)
ax2.legend(loc = 'upper left')
ax1.legend(loc = 'upper right')
plt.show()


#import kendalltau to be able to find correlation 
from scipy.stats import kendalltau

#create an empty array of correlation coefficients that is the size of the number of files
corr_coeffs = [0] * np.size(filenames)

#create an empty array of pvalues that is the size of the number of files
pvalues = [0] * np.size(filenames)

#define count to loop through the arrays defined above 
count = 0

#write a loop to find the kendall tau correlation coefficient 
for name in filenames:
    #similar to before, we get the key by removing the .csv ending from the filename
    key = name.split('.')[0]
    #extract just nitrogen data
    Nitrogen = nutrient_data[key]['Nitrate_Nitrite']
    #extract just TDP data
    Phosphorus = nutrient_data[key]['Total_Dissolved_Phosphorus']
    #calculate the correlation
    corr_coeff, pvalue = kendalltau(Nitrogen, Phosphorus)
    print(f"The Kendall's \N{GREEK SMALL LETTER TAU} coefficient is {corr_coeff:.4f} and the pvalue is {pvalue:.5f}" 
          f"at Site {key}.")
    #add to an array of correlation coefficients for printing purposes later
    corr_coeffs[count] = corr_coeff
    #do the same for pvalues
    pvalues[count] = pvalue
    #increment count
    count = count + 1

The Kendall's τ coefficient is 0.2163 and the pvalue is 0.00026at Site WE2.
The Kendall's τ coefficient is 0.0298 and the pvalue is 0.61924at Site WE4.
The Kendall's τ coefficient is 0.4376 and the pvalue is 0.00000at Site WE6.
The Kendall's τ coefficient is 0.2457 and the pvalue is 0.00003at Site WE8.
The Kendall's τ coefficient is 0.3694 and the pvalue is 0.00004at Site WE9.
The Kendall's τ coefficient is 0.1998 and the pvalue is 0.00378at Site WE12.
The Kendall's τ coefficient is -0.0106 and the pvalue is 0.88122at Site WE13.


#extract data for just 2015 for site WE12 using the dates datetime column
WE12_2015 = [row for row in nutrient_data['WE12'] if row['\ufeffDate'].year == 2015]

#define a function that will allow us to extract nutrient data from the data subsets
def extract_values(dataset, nutrient):
    '''
    Allows us to easily extract a specific subset of data based on the column name, in the parameters 
    this is 'nutrient'. Returns all data in a specified dataset for that nutrient
    '''
    return([row[nutrient] for row in dataset])

#Extract just the dates from the 2015 sub-data (this will be the x-axis)
dates_2015 = extract_values( WE12_2015, '\ufeffDate')
    
#Extract just nutrient info from the 2015 sub-data (this will be the y-axis)
Nvalues_2015 = extract_values(WE12_2015,'Nitrate_Nitrite')

#define a figure 
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize = (10, 12))

#plot the nitrogen data as a line with dots
ax1.plot(dates_2015, Nvalues_2015, 'o-', color = 'forestgreen')

#as for 2015, extract data for just 2016 for site WE12 using the dates datetime column 
WE12_2016 = [row for row in nutrient_data['WE12'] if row['\ufeffDate'].year == 2016]

#Extract just the dates from the 2016 sub-data (this will be the x-axis)
dates_2016 = extract_values( WE12_2016, '\ufeffDate')

#Extract just nutrient info from the 2016 sub-data (this will be the y-axis)
Nvalues_2016 = extract_values(WE12_2016,'Nitrate_Nitrite')

#plot the nitrogen data as a line with dots
ax2.plot(dates_2016, Nvalues_2016, 'o-', color = 'forestgreen')

#add plot elements to ax1 (labels, legend, title) 
ax1.set_title('Nitrate/Nitrite, 2015', fontsize='x-large')
ax1.set_xlabel('Year', fontsize='x-large')
ax1.set_ylabel('N/N Concentrations [$mg N/L$]')
ax1.tick_params(axis = 'x', rotation = 45)
ax1.xaxis.set_major_formatter(mdates.DateFormatter("%m-%d"))

#add plot elements to ax2 (labels, legend, title) 
ax2.set_title('Nitrate/Nitrite, 2016', fontsize='x-large')
ax2.set_xlabel('Year', fontsize='x-large')
ax2.set_ylabel('N/N Concentrations [$mg N/L$]')
ax2.tick_params(axis = 'x', rotation=45)
ax2.xaxis.set_major_formatter(mdates.DateFormatter("%m-%d"))

#Extract just TDP info from the 2015 sub-data (this will be the y-axis)
Pvalues_2015 = extract_values(WE12_2015,'Total_Dissolved_Phosphorus')

#plot the TDP data as lines with dots 
ax3.plot(dates_2015, Pvalues_2015, 'o-', color = 'turquoise')

#Extract just TDP info from the 2016 sub-data (this will be the y-axis)
Pvalues_2016 = extract_values(WE12_2016, 'Total_Dissolved_Phosphorus')

#plot the TDP data as lines with dots 
ax4.plot(dates_2016, Pvalues_2016, 'o-', color = 'turquoise')

#add plot elements to ax3 (labels, legend, title) 
ax3.set_title('Total Dissolved Phosphorus, 2015', fontsize='x-large')
ax3.set_xlabel('Year', fontsize='x-large')
ax3.set_ylabel('TDP concentrations [$\mu$g P/L]')
ax3.tick_params(axis = 'x', rotation = 45)
ax3.xaxis.set_major_formatter(mdates.DateFormatter("%m-%d"))

#add plot elements to ax4 (labels, legend, title) 
ax4.set_title('Total Dissolved Phosphorus, 2016', fontsize='x-large')
ax4.set_xlabel('Year', fontsize='x-large')
ax4.set_ylabel('TDP Concentrations [$\mu$g P/L]')
ax4.tick_params(axis = 'x', rotation = 45)
ax4.xaxis.set_major_formatter(mdates.DateFormatter("%m-%d"))

#give the figure an overarching title 
fig.suptitle('Nutrient Concentrations at Site WE12 During Significant (2015) versus Mild (2016) HABs Years', fontsize = 16)

#adjust the figure layout
plt.tight_layout(h_pad = 1.5, w_pad = 1.5)


#statistics for the above plots 

#define an array of the data plotted above 
data = [Nvalues_2015, Nvalues_2016,  Pvalues_2015, Pvalues_2016]

#define names that can be used in print statements for above data
names = ["Nitrogen 2015 (mg/L)","Nitrogen 2016 (mg/L)", "TDP 2015 (\u00B5g P/L)", "TDP 2016 (\u00B5g P/L)"]

#counter
i = 0

#loop through the data and find the max, min, and average value
#print out these values
for x in data: 
    max = np.max(x)
    min = np.min(x)
    average = np.average(x)
    print(f'The maximum value in {names[i]} is {np.round(max, 3)}. The minimum '
           f'value is {np.round(min, 3)}. The average value is {np.round(average, 3)}.')
    i = i + 1 
    
#add a space
print()

#print the Kendall Tau correlation and pvalue for 2015 and 2016 nutrients
corr_coeff_2015, pvalue_2015 = kendalltau(Nvalues_2015, Pvalues_2015)

corr_coeff_2016, pvalue_2016 = kendalltau(Nvalues_2016, Pvalues_2016)

print(f'The Kendall \N{GREEK SMALL LETTER TAU} correlation between nitrogen and phosphorus in the 2015 '
      f'(significant) HABs year was {np.round(corr_coeff_2015, 3)}.'
f'with a pvalue of {np.round(pvalue_2015, 5)}.')
print(f'The Kendall \N{GREEK SMALL LETTER TAU} correlation between nitrogen and phosphorus in the 2016 '
      f'(mild) HABs year was {np.round(corr_coeff_2016, 3)} with a pvalue of {np.round(pvalue_2016, 3)}.')

The maximum value in Nitrogen 2015 (mg/L) is 6.25. The minimum value is 0.002. The average value is 1.646.
The maximum value in Nitrogen 2016 (mg/L) is 1.26. The minimum value is 0.0. The average value is 0.279.
The maximum value in TDP 2015 (µg P/L) is 98.27. The minimum value is 4.92. The average value is 44.53.
The maximum value in TDP 2016 (µg P/L) is 32.0. The minimum value is 3.9. The average value is 12.595.

The Kendall τ correlation between nitrogen and phosphorus in the 2015 (significant) HABs year was 0.67.with a pvalue of 4e-05.
The Kendall τ correlation between nitrogen and phosphorus in the 2016 (mild) HABs year was 0.048 with a pvalue of 0.778.


#WE6 N/N and TDP concentrations in 2015 (HAB year) vs 2012 (non-HAB year)
WE6_2015 = [row for row in nutrient_data['WE6'] if row[0].year == 2015]

#Extract datetime values and measurements
dates6_2015 = np.array([row[0] for row in WE6_2015])
N6values_2015 = [row['Nitrate_Nitrite'] for row in WE6_2015]

fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize = (10, 12))
ax1.scatter(dates6_2015, N6values_2015,  color = 'forestgreen')
ax1.plot(dates6_2015, N6values_2015,  color = 'forestgreen')

WE6_2012 = [row for row in nutrient_data['WE6'] if row[0].year == 2012]
# Extract datetime values and measurements
dates6_2012 = np.array([row[0] for row in WE6_2012])
N6values_2012 = [row['Nitrate_Nitrite'] for row in WE6_2012]

ax2.scatter(dates6_2012, N6values_2012, color = 'forestgreen')
ax2.plot(dates6_2012, N6values_2012, color = 'forestgreen')

ax1.set_title('Nitrate/Nitrite, 2015', fontsize='x-large')
ax1.set_xlabel('Year', fontsize='x-large')
ax1.set_ylabel('N/N Concentrations [$mg N/L$]', fontsize='large')
ax1.tick_params(axis = 'x', rotation = 45)
ax1.xaxis.set_major_formatter(mdates.DateFormatter("%m-%d"))
ax2.set_title('Nitrate/Nitrite, 2012', fontsize='x-large')
ax2.set_xlabel('Year', fontsize='x-large')
ax2.set_ylabel('N/N Concentrations[$mg N/L$]', fontsize='large')
ax2.tick_params(axis = 'x', rotation=45)
ax2.xaxis.set_major_formatter(mdates.DateFormatter("%m-%d"))

P6values_2015 = [row['Total_Dissolved_Phosphorus'] for row in WE6_2015]
ax3.scatter(dates6_2015, P6values_2015, color = 'turquoise')
ax3.plot(dates6_2015, P6values_2015, color = 'turquoise')
ax3.tick_params(axis = 'x', rotation = 45)

P6values_2012 = [row['Total_Dissolved_Phosphorus'] for row in WE6_2012]
ax4.scatter(dates6_2012, P6values_2012, color = 'turquoise')
ax4.plot(dates6_2012, P6values_2012, color = 'turquoise')
ax4.tick_params(axis = 'x', rotation = 45)

ax3.set_title('Total Dissolved Phosphorus, 2015', fontsize='x-large')
ax3.set_xlabel('Year', fontsize='x-large')
ax3.set_ylabel('TDP Concentrations [$\mu$g P/L]', fontsize='large')
ax3.xaxis.set_major_formatter(mdates.DateFormatter("%m-%d"))
ax4.set_title('Total Dissolved Phosphorus, 2016', fontsize='x-large')
ax4.set_xlabel('Year', fontsize='x-large')
ax4.set_ylabel('TDP Concentrations [$\mu$g P/L]', fontsize='large')
ax4.xaxis.set_major_formatter(mdates.DateFormatter("%m-%d"))

fig.suptitle('Nutrient Concentrations at Site WE6 During Significant (2015) versus Mild (2012) HABs Event Years', fontsize = 15)
plt.tight_layout(h_pad = 1.5, w_pad = 1.5)


#WE6 N/N and TDP concentrations in 2015 (HAB year) vs 2016 (non-HAB year)
WE6_2015 = [row for row in nutrient_data['WE6'] if row[0].year == 2015]

#Extract datetime values and measurements
dates6_2015 = np.array([row[0] for row in WE6_2015])
N6values_2015 = [row['Nitrate_Nitrite'] for row in WE6_2015]

fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize = (10, 12))
ax1.scatter(dates6_2015, N6values_2015,  color = 'forestgreen')
ax1.plot(dates6_2015, N6values_2015,  color = 'forestgreen')

WE6_2016 = [row for row in nutrient_data['WE6'] if row[0].year == 2016]
# Extract datetime values and measurements
dates6_2016 = np.array([row[0] for row in WE6_2016])
N6values_2016 = [row['Nitrate_Nitrite'] for row in WE6_2016]

ax2.scatter(dates6_2016, N6values_2016, color = 'forestgreen')
ax2.plot(dates6_2016, N6values_2016, color = 'forestgreen')

ax1.set_title('Nitrate/Nitrite Concentrations in 2015')
ax1.set_xlabel('Dates [$year$]')
ax1.set_ylabel('N/N [$mg N/L$]')
ax1.tick_params(axis = 'x', rotation = 45)
ax2.set_title('Nitrate/Nitrite Concentrations in 2016)')
ax2.set_xlabel('Dates [$year$]')
ax2.set_ylabel('N/N [$mg N/L$]')
ax2.tick_params(axis = 'x', rotation=45)

P6values_2015 = [row['Total_Dissolved_Phosphorus'] for row in WE6_2015]
ax3.scatter(dates6_2015, P6values_2015, color = 'turquoise')
ax3.plot(dates6_2015, P6values_2015, color = 'turquoise')
ax3.tick_params(axis = 'x', rotation = 45)

P6values_2016 = [row['Total_Dissolved_Phosphorus'] for row in WE6_2016]
ax4.scatter(dates6_2016, P6values_2016, color = 'turquoise')
ax4.plot(dates6_2016, P6values_2016, color = 'turquoise')
ax4.tick_params(axis = 'x', rotation = 45)

ax3.set_title('Total Dissolved Phosphorus Concentrations in 2015')
ax3.set_xlabel('Dates [$year$]')
ax3.set_ylabel('TDP [$\mu$g P/L]')
ax4.set_title('Total Dissolved Phosphorus Concentrations in 2016')
ax4.set_xlabel('Dates [$year$]')
ax4.set_ylabel('TDP [$\mu$g P/L]')

fig.suptitle('Nutrient Concentrations at Site WE6 in Western Lake Erie Basin During Significant (2015) versus Mild (2016) HABs Years', fontsize = 15)
plt.tight_layout(h_pad = 1.5, w_pad = 1.5)


#we are first going to compare WE6 to WE13, WE6 is the larger dataset

#upload the reformatted WE6 file so all data points match WE13 file
WE6_reformat = nutrient_data_file_reader('WE6_Reformat.csv')


#extract the nitrogen data
N_WE6 = WE6_reformat['Nitrate_Nitrite']

#extract the ammonia data 
A_WE6 = WE6_reformat['Ammonia']

#extract the TDP data
TDP_WE6 = WE6_reformat['Total_Dissolved_Phosphorus']

#extract the TP data
P_WE6 = WE6_reformat['TotalPhosphorus']

#extract the nitrogen data
N_WE13 = nutrient_data['WE13']['Nitrate_Nitrite']

#extract the ammonia data 
A_WE13 = nutrient_data['WE13']['Ammonia']

#extract the TDP data
TDP_WE13 = nutrient_data['WE13']['Total_Dissolved_Phosphorus']

#extract the TP data
P_WE13 = nutrient_data['WE13']['TotalPhosphorus']

#stack all the results together into one result set to be used for comparison
results = np.column_stack((N_WE6, A_WE6 , TDP_WE6, P_WE6, N_WE13,A_WE13, TDP_WE13, P_WE13))

#transpose the column stack result 
dataT = results.T

# Get the correlation coeff matrix:
#we expect the diagonal to be 1s bc it is the correlation of a variable with itself
cc = np.corrcoef(dataT)

# Create a figure and add an axes object:
fig, (ax1, ax2) = plt.subplots(2, 1, figsize = (5, 10))

#adjust spacing between subplots
plt.subplots_adjust(hspace=0.8)

#Create the "pcolor" object while plotting it to the axes.
#Use our pcolor object to create a colorbar.
pcol = ax1.pcolor(cc, cmap='seismic', vmin=-1., vmax=1.)
cbar = plt.colorbar(pcol, ax=ax1)
cbar.set_label('Correlation')

#Set labels and location of labels

#N&N = Nitrate/Nitrite
#A = Ammonia
#TDP = total dissolved phosphorus
#TP = total phosphorus
labels = ['N&N (WE6)', 'A (WE6)', 'TDP (WE6)','TP (WE6)', 'N&N (WE13)', 'A (WE13)', 'TDP (WE13)', 'TP (WE13)']
locs = np.arange(0, 8, 1.) + .5

#Make everything line up properly
ax1.set_xticks(locs, labels, rotation= 'vertical')
ax1.set_yticks(locs, labels)

#Title
ax1.set_title('Nutrient Correlation for Limnetic Zone vs Littoral Zone Sites: Spatially Separated', pad=14)
#A Limnetic Zone is towards the middle of a lake, far into the deeper part of a body of water --> further from land!

##NEXT CODE BLOCK IS ALMOST IDENTICAL TO ABOVE FOR DIFFERENT SITES (WE6 and WE9)

#next we are going to compare WE6 to WE9, WE6 is the larger dataset
#upload the reformatted WE6 file so all data points match WE9 file
WE6_reformat_2 = np.genfromtxt('WE6_Reformat_2.csv', delimiter=',', converters={0: tconvert}, 
                               encoding='utf-8', names=True, dtype=None)

#extract the nitrogen data
N_WE6 = WE6_reformat_2['Nitrate_Nitrite']

#extract the ammonia data 
A_WE6 = WE6_reformat_2['Ammonia']

#extract the TDP data
TDP_WE6 = WE6_reformat_2['Total_Dissolved_Phosphorus']

#extract the TP data
P_WE6 = WE6_reformat_2['TotalPhosphorus']

#Do the same extractions for WE()
N_WE9 = nutrient_data['WE9']['Nitrate_Nitrite']
A_WE9 = nutrient_data['WE9']['Ammonia']
TDP_WE9 = nutrient_data['WE9']['Total_Dissolved_Phosphorus']
P_WE9 = nutrient_data['WE9']['TotalPhosphorus']

#add the results together using column stack
results = np.column_stack((N_WE6, A_WE6 , TDP_WE6, P_WE6, N_WE9,A_WE9, TDP_WE9, P_WE9))

#transpose the results 
dataT = results.T

#Get the correlation coeff matrix:
#we expect the diagonal to be 1s bc it is the correlation of a variable with itself
cc = np.corrcoef(dataT)


#Create the "pcolor" object while plotting it to the axes.
#Use our pcolor object to create a colorbar.
pcol = ax2.pcolor(cc, cmap = 'seismic', vmin = -1., vmax = 1.)
cbar = plt.colorbar(pcol, ax = ax2)
cbar.set_label('Correlation')

#Set labels and location
labels = ['N&N (WE6)', 'A (WE6)', 'TDP (WE6)','TP (WE6)', 'N&N (WE9)', 'A (WE9)', 'TDP (WE9)', 'TP (WE9)']
locs = np.arange(0, 8, 1.) + .5

#Make sure everything lines up 
ax2.set_xticks(locs, labels, rotation = 'vertical')
ax2.set_yticks(locs, labels)

#Title our axes:
ax2.set_title('Nutrient Data Correlation for Littoral Zone Sites: Spatially Similar', pad=14)
#A Littoral Zone is closest to the coast of a body of water, so near land!

plt.show()


#these data files were exactly the same size with matching dates, so we did not have to do any manipulation here

#extract nitrate/nitrite
n_N_WE6 = nutrient_data['WE6']['Nitrate_Nitrite']

#extract ammonia
n_A_WE6 = nutrient_data['WE6']['Ammonia']

#extract TDP
n_TDP_WE6 = nutrient_data['WE6']['Total_Dissolved_Phosphorus']

#extract TP
n_P_WE6 = nutrient_data['WE6']['TotalPhosphorus']

#Do the same nutrient extraction for site WE8
N_WE8 = nutrient_data['WE8']['Nitrate_Nitrite']
A_WE8 = nutrient_data['WE8']['Ammonia']
TDP_WE8 = nutrient_data['WE8']['Total_Dissolved_Phosphorus']
P_WE8 = nutrient_data['WE8']['TotalPhosphorus']


#add all the results together into one contiguous result set 
results = np.column_stack((n_N_WE6, n_A_WE6 , n_TDP_WE6, n_P_WE6, N_WE8,A_WE8, TDP_WE8, P_WE8))

#transpose the data
dataT = results.T

# Get the correlation coeff matrix:
#we expect the diagonal to be 1s bc it is the correlation of a variable with itself
cc = np.corrcoef(dataT)

#define a figure with 1 subplot
fig, (ax2) = plt.subplots( 1, figsize = (7, 5))

# Create the "pcolor" object while plotting it to the axes.
# Use our pcolor object to create a colorbar.
pcol = ax2.pcolor(cc, cmap='seismic', vmin=-1., vmax=1.)
cbar = plt.colorbar(pcol, ax=ax2)
cbar.set_label('Correlation')

# Set labels and label locations
labels = ['N/N (WE6)', 'A (WE6)', 'TDP (WE6)','TP (WE6)', 'N/N (WE8)', 'A (WE8)', 'TDP (WE8)','TP (WE8)']
locs = np.arange(0, 8, 1.) + .5

# Now, set our ticks and labels
ax2.set_xticks(locs, labels, rotation= 'vertical')
ax2.set_yticks(locs, labels)

# Title our axes:
ax2.set_title('Nutrient Correlation for Littoral Zone Sites: Spatially Separated', pad=15)
plt.show()


#define the variable names 
varnames = ['Date', 'WaterTemps']

#Load water data without converters
lake_erie_temps = np.genfromtxt('Water_Temps_2012_2018.csv', delimiter = ',', names = varnames, dtype = None, 
                                encoding = 'utf-8', skip_header = 1)

#convert to datetime objects 
temp_dates = [dt.datetime.strptime(date_str, '%m/%d/%y') for date_str in lake_erie_temps['Date']]

#create an object to hold all water temperatures
water_temps  = lake_erie_temps['WaterTemps']

#create a figure
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 6))

#adjust spacing
plt.subplots_adjust(hspace=0.5)

#plot all water temp data 
ax1.plot(temp_dates, water_temps)

#add labels and title
ax1.set_ylabel('Water Temp (ºC)')
ax1.set_xlabel('Year')
ax1.set_title('Daily Water Temperature in Lake Erie 2012-2018')

#define the number of points in the data set
N = water_temps.size

#find the FFT of the water data
data_freq = fft(water_temps)

#create a frequency array for the data set in cycles per year
freqs_days = fftfreq(N, 1/365)

#create an amplitude array 
amps = np.abs(data_freq)

# Normalize and plot the positive data 
ax2.plot(freqs_days[1:N // 2], amps[1:N // 2]/N)

# Add labels and title
ax2.set_xlabel('Frequency [$cycles/year$]')
ax2.set_ylabel("Amplitude [ºC]")
ax2.set_title("Amplitude Spectrum of Water Temperatures in Lake Erie 2012-2018")
ax2.set_xlim(0, 5)

# Find the index of the highest frequency 
max_index = np.argmax(amps[1:N // 2]/N)

# Find the frequency at that index
fmax = freqs_days[max_index]

# Use the max frequency to find the max amplitude 
amp_highest_freq = amps[max_index]

# Find the period associated with fmax
max_period = (1/fmax)

print(f'The highest frequency in this dataset is {fmax} and the associted period is '
      f'{max_period}.')

The highest frequency in this dataset is 0.8744010951403148 and the associted period is 1.1436399217221136.


fig, (ax1) = plt.subplots(1, figsize=(10, 5))
#maybe plot each site color the same as the map?
count = 0
colors = ['maroon', 'orangered', 'orange', 'gold', 'green', 'blue', 'blueviolet', 'magenta']
for name in filenames:
    key = name.split('.')[0]
    color = colors[count]
    count = count + 1
    ax1.plot(nutrient_data[key]['\ufeffDate'], nutrient_data[key]['Nitrate_Nitrite'], label = key, color = color)


ax1.legend(loc='upper left', bbox_to_anchor=(-0.2, 1))
ax1.set_xlabel('Year', fontsize='x-large')
ax1.set_ylabel('N/N Concentrations [$mg N/L$]', fontsize='large')
ax1.set_title('Nitrate and Nitrite Concentrations at Sites in Western Lake Erie Basin 2012-2018', fontsize = 15)

ax2 = ax1.twinx()
ax2.set_ylabel('Water Temp (ºC)')
ax2.plot(temp_dates[360:], water_temps[360:], color = 'fuchsia', label = 'Temps', linewidth= 3, alpha = 0.8)
ax2.legend(loc='upper left')
fig.tight_layout()


#create a figure with two subplots
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 5))

#array of colors we want our lines plotted as 
colors = ['maroon', 'orangered', 'orange', 'gold', 'green', 'blue', 'blueviolet', 'magenta']
#define a variable to help iterature through the colors array 
count = 0

#same loop as previous code block, but we are limiting to 2016 data 
for name in filenames:
    key = name.split('.')[0]
    color = colors[count]
    count += 1
    
    #FILTERING
    #extract all dates at a specific site
    dates = nutrient_data[key]['\ufeffDate']
    
    #extract all N/N data at a specific site 
    N_concentrations = nutrient_data[key]['Nitrate_Nitrite']
    
    #extract all P data at a specific site 
    P_concentrations = nutrient_data[key]['Total_Dissolved_Phosphorus']
    
    #only get the indices of dates that are within the desired year (2106)
    year_2016_indices = [i for i, date in enumerate(dates) if date.year == 2016]
    
    #use the indices we just found to go back through the date array and include only those dates
    dates_2016 = [dates[i] for i in year_2016_indices]
    
    #only extract nutrient concentrations for the desired year (2016)
    N_concentrations_2016 = [N_concentrations[i] for i in year_2016_indices]
    P_concentrations_2016 = [P_concentrations[i] for i in year_2016_indices]
    
    #plot the sub-data set 
    ax1.plot(dates_2016, N_concentrations_2016 , label=key, color=color)
    ax2.plot(dates_2016, P_concentrations_2016, label=key, color=color)

#add plot elements - legend, title, labels
ax1.legend(loc='upper left', bbox_to_anchor=(-0.2, 1))
ax1.set_xlabel('Year', fontsize='x-large')
ax1.set_ylabel('N/N Concentrations [$mg N/L$]', fontsize='large')
ax1.set_title('Nitrate and Nitrite Concentrations Compared to Tempeerature at All Sites in 2016', fontsize=15)

ax2.legend(loc='upper left', bbox_to_anchor=(-0.2, 1))
ax2.set_xlabel('Year', fontsize='x-large')
ax2.set_ylabel('TDP Concentrations [$\mu$g P/L]', fontsize='large')
ax2.set_title('TDP Concentrations Compared to Tempeerature at All Sites in 2016', fontsize=15)


#extract the dates from the temperature dataset that are within the time period we are interested in 
#(5/1/2016 - 12/31/2016)
temp_dates_2015 = [date for date in temp_dates if dt.datetime(2016, 5, 1) <= date <= dt.datetime(2016, 12, 31)]

#do the same for the actual temperature data
water_temps_2015 = [water_temps[i] for i, date in enumerate(temp_dates) if dt.datetime(2016, 5, 1) <= date
                    <= dt.datetime(2016, 12, 31)]

#create a twin axes for the first subplot
ax3 = ax1.twinx()
#Plot water temperatures for those 2016 dates on first subplot
ax3.plot(temp_dates_2015, water_temps_2015, color='fuchsia', label='Temps', linewidth=3)

#create a twin axes for the second subplot
ax4 = ax2.twinx()

#Plot water temperatures for those 2016 dates on second subplot
ax4.plot(temp_dates_2015, water_temps_2015, color='fuchsia', label='Temps', linewidth=3)

#add labels
ax3.set_ylabel('Water Temp (ºC)')
ax3.legend(loc='upper left')
ax4.set_ylabel('Water Temp (ºC)')
ax4.legend(loc='upper left')

#adjust layout
fig.tight_layout()
plt.show()


# import runoff data:
# convert first 2 columns to datetime object
tconvert = lambda x: dt.datetime.strptime(str(x), '%Y,  %m,')

# create list of new variable names
varnames = ['datetime', 'flow']

# use genfromtxt to read in data
# set variable names, skip 3 lines for header, extract only relevant columns
# specify column widths & set converter to extract first 2 columns into datetime object!
runoff_data = np.genfromtxt('LakeErie-runoff.csv', names=varnames, delimiter=(9,7), skip_header=3, usecols=(0,1),
                     encoding='utf-8', converters={0:tconvert}, dtype=None)

# data clean-up: delete days where flow value = -999.9
runoff_data = runoff_data[runoff_data['flow'] != -999.9]


# import wastewater discharge data:
# convert first column to datetime object
tconvert = lambda x: dt.datetime.strptime(str(x), '%m/%d/%Y')

# create list of new variable names
varnames = ['datetime', 'final_ammonia', 'secondary_TP', 'secondary_NH3']

# use genfromtxt to read in data
# set converter for the zero-th column to create datetime object, set variable names, skip header
ww_data = np.genfromtxt('GLWA Effluent Nutrients 2012-2018.csv', names=varnames, delimiter=',', skip_header=1,
                     encoding='utf-8', converters={0:tconvert}, dtype=None)

# some 'secondary_TP' measurements are below threshold; set to 0!
ww_data['secondary_TP'][ww_data['secondary_TP'] == "<0.02"] = 0
ww_data['secondary_TP'] = ww_data['secondary_TP'].astype(float)


# filter data! only keep values for which dates match. (i.e., nutrient data and runoff data should be the same length!)

# define a function to filter data down to size
def time_filter(dataset, site):
    '''
    This function filters a dataset by date, keeping ONLY rows with a datetime present in the specified
    nutrient dataset. Returns a filtered data array of the SAME SIZE as the nutrient data for the particular
    site, with MATCHING, ONE-TO-ONE datetimes. This allows for scatter plotting and bivariate analysis between
    both data arrays.
    
    ========
    args:
    data = longer numpy array containing data from extra datetimes; this dataset is to be filtered/cut to size
    site = name of desired nutrient data site; the supplied longer data array will be cut against the data for the site
    
    ========
    example usage:
    data_filtered = timefilter(data, 'WE6')
    
    # resulting data_filtered is now the exact same size as nutrient_data['WE6'], with the exact same 'datetime' values!
    # they can now easily be plotted against one another / compared
    
    '''
    
    # initialize array for boolean index
    index = []
    
    # iterate through datetimes of longer dataset, seeing if datetime matches with nutrient measurement datetime
    # store (True/False) value in boolean array
    for dtime in dataset['datetime']:
        index.append(dtime in nutrient_data[site]['\ufeffDate'])
    
    # index data by resulting boolean array. return!
    return dataset[index]


# designate site
site = 'WE4'

# filter out wastewater data based on datetimes of given site
ww_data_filtered = time_filter(ww_data, site)

# create figure, axes objects
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(12,10))
# adjust vertical spacing
plt.subplots_adjust(hspace=0.4)

# plot/compare a different nutrient on each subplot!
ax1.scatter(nutrient_data[site]['TotalPhosphorus']/1000, ww_data_filtered['final_ammonia'], color='rebeccapurple')
ax1.set_title(f"Wastewater Ammonia vs. TP at site {site}", fontsize='x-large')
ax1.set_xlabel(f"TP Concentrations, {site} [$mg/L$]", fontsize='large')
ax1.set_ylabel("Effluent Ammonia [$mg/L$]", fontsize='large')

ax2.scatter(nutrient_data[site]['Total_Dissolved_Phosphorus']/1000, ww_data_filtered['final_ammonia'], color='rebeccapurple')
ax2.set_title(f"Wastewater Ammonia vs. TDP at site {site}", fontsize='x-large')
ax2.set_xlabel(f"TDP Concentrations, {site} [$mg/L$]", fontsize='large')
ax2.set_ylabel("Effluent Ammonia [$mg/L$]", fontsize='large')

ax3.scatter(nutrient_data[site]['Ammonia']/1000, ww_data_filtered['final_ammonia'], color='rebeccapurple')
ax3.set_title(f"Wastewater Ammonia vs. Ammonia at site {site}", fontsize='x-large')
ax3.set_xlabel(f"Ammonia Concentrations, {site}", fontsize='large')
ax3.set_ylabel("Effluent Ammonia [$mg/L$]", fontsize='large')

ax4.scatter(nutrient_data[site]['Nitrate_Nitrite'], ww_data_filtered['final_ammonia'], color='rebeccapurple')
ax4.set_title(f"Wastewater Ammonia vs. N/N at site {site}", fontsize='x-large')
ax4.set_xlabel(f"N/N Concentrations, {site} [$mg/L$]", fontsize='large')
ax4.set_ylabel("Effluent Ammonia [$mg/L$]", fontsize='large')
plt.show()

# iterate through plots, calculate Kendall's tau for each
print("Kendall's \N{GREEK SMALL LETTER TAU} coefficient:")
for meas in ['TotalPhosphorus', 'Total_Dissolved_Phosphorus', 'Ammonia', 'Nitrate_Nitrite']:
    tau, p_value = kendalltau(nutrient_data[site][meas], ww_data_filtered['final_ammonia'])
    print(f"Wastewater ammonia vs. {meas} = {tau:.4f}, p = {p_value:.4f}")

Kendall's τ coefficient:
Wastewater ammonia vs. TotalPhosphorus = -0.1285, p = 0.0317
Wastewater ammonia vs. Total_Dissolved_Phosphorus = -0.0964, p = 0.1072
Wastewater ammonia vs. Ammonia = -0.0849, p = 0.1558
Wastewater ammonia vs. Nitrate_Nitrite = 0.0793, p = 0.1875


# take monthly averages for nutrient measurements 
def nutrient_mean(site):
    '''
    This function returns an array of MONTHLY AVERAGES for nutrient data, from the original nutrient 
    concentration data at the specified site. Datetimes for monthly averages are set to the first of
    the month.
    '''
    
    avg = [dt.datetime(1, 1, 1), 0, 0, 0, 0]
    for year in range(2012, 2019, 1):
        for month in range(1,13):
            temp = nutrient_data[site][0]
            for i, dtime in enumerate(nutrient_data[site]['\ufeffDate']):
                if (dtime.month == month) & (dtime.year == year):
                    temp = np.vstack((temp, nutrient_data[site][i]))
            if(temp.size > 1):
                temp = temp[1:]       
                month_row = np.array([dt.datetime(year, month, 1), temp['TotalPhosphorus'].mean(), 
                                      temp['Total_Dissolved_Phosphorus'].mean(), temp['Ammonia'].mean(), 
                                      temp['Nitrate_Nitrite'].mean()])
                avg = np.vstack((avg, month_row))

    avg = avg[1:]
    return avg


# for site WE2:
site = 'WE2'

# obtain monthly averages of relevant nutrient data
nutrient_monthly = nutrient_mean(site)

# shorten runoff data to range (use same algorithm as in time_filter() function)
index = []
for dtime in runoff_data['datetime']:
    index.append(dtime in nutrient_monthly[:,0])
runoff_shortened = runoff_data[index]

# plot
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(12,9))
plt.subplots_adjust(hspace=0.4)

ax1.scatter(nutrient_monthly[:,1], runoff_shortened['flow'], color="rebeccapurple")
ax1.set_title(f"Runoff vs. TP at site {site}", fontsize='x-large')
ax1.set_xlabel(f"TP Concentrations, {site} [$μg/L$]", fontsize='large')
ax1.set_ylabel("Monthly Runoff [$m^3/d$]", fontsize='large')
# find fit coefficients
coefs1 = np.polyfit(np.array(nutrient_monthly[:,1], dtype=float), runoff_shortened['flow'], 1)
# create a new x-array, obtain polyfit y-values
y1_fit = np.polyval(coefs1, nutrient_monthly[:,1])
ax1.plot(nutrient_monthly[:,1], y1_fit, color='rebeccapurple')

ax2.scatter(nutrient_monthly[:,2], runoff_shortened['flow'], color="rebeccapurple")
ax2.set_title(f"Runoff vs. TDP at site {site}", fontsize='x-large')
ax2.set_xlabel(f"TDP Concentrations, {site} [$μ/L$]", fontsize='large')
ax2.set_ylabel("Monthly Runoff [$m^3/d$]", fontsize='large')
coefs2 = np.polyfit(np.array(nutrient_monthly[:,2], dtype=float), runoff_shortened['flow'], 1)
# create a new x-array, obtain polyfit y-values
y2_fit = np.polyval(coefs2, nutrient_monthly[:,2])
ax2.plot(nutrient_monthly[:,2], y2_fit, color='rebeccapurple')

ax3.scatter(nutrient_monthly[:,3], runoff_shortened['flow'], color="rebeccapurple")
ax3.set_title(f"Runoff vs. Ammonia at site {site}", fontsize='x-large')
ax3.set_xlabel(f"Ammonia Concentrations, {site} [$μ/L$]", fontsize='large')
ax3.set_ylabel("Monthly Runoff [$m^3/d$]", fontsize='large')
coefs3 = np.polyfit(np.array(nutrient_monthly[:,3], dtype=float), runoff_shortened['flow'], 1)
# create a new x-array, obtain polyfit y-values
y3_fit = np.polyval(coefs1, nutrient_monthly[:,3])
ax3.plot(nutrient_monthly[:,3], y3_fit, color='rebeccapurple')

ax4.scatter(nutrient_monthly[:,4], runoff_shortened['flow'], color="rebeccapurple")
ax4.set_title(f"Runoff vs. N/N at site {site}", fontsize='x-large')
ax4.set_xlabel(f"N/N Concentrations, {site} [$mg/L$]", fontsize='large')
ax4.set_ylabel("Monthly Runoff [$m^3/d$]", fontsize='large')
coefs4 = np.polyfit(np.array(nutrient_monthly[:,4], dtype=float), runoff_shortened['flow'], 1)
# create a new x-array, obtain polyfit y-values
y4_fit = np.polyval(coefs4, nutrient_monthly[:,4])
ax4.plot(nutrient_monthly[:,4], y4_fit, color='rebeccapurple')

plt.show()
    
print("Kendall \N{GREEK SMALL LETTER TAU} coefficient:")
for i, meas in enumerate(['TotalPhosphorus', 'Total_Dissolved_Phosphorus', 'Ammonia', 'Nitrate_Nitrite']):
    tau, p_value = kendalltau(nutrient_monthly[:,i+1].astype(float), runoff_shortened['flow'])
    print(f"Runoff vs. {meas} at {site} = {tau:.4f}, p = {p_value:.4f}")

Kendall τ coefficient:
Runoff vs. TotalPhosphorus at WE2 = 0.2558, p = 0.0156
Runoff vs. Total_Dissolved_Phosphorus at WE2 = 0.3533, p = 0.0008
Runoff vs. Ammonia at WE2 = 0.2447, p = 0.0207
Runoff vs. Nitrate_Nitrite at WE2 = 0.5216, p = 0.0000


# create a matrix with all correlation coef.s and Kendall tau coef.s (including all sites vs. all measured variables)
CC_all = np.zeros((4, 7))
Kendall_all = np.zeros((4, 7))

# iterating through sites
for j, site in enumerate(['WE2', 'WE4', 'WE6', 'WE8', 'WE9', 'WE12', 'WE13']):
    nutrient_monthly = nutrient_mean(site)

    # shorten runoff data to range (use same algorithm as in earlier function)
    index = []
    for dtime in runoff_data['datetime']:
        index.append(dtime in nutrient_monthly[:,0])
    runoff_shortened = runoff_data[index]

    # iterate through each nutrient metric in each site
    for i, meas in enumerate(['TotalPhosphorus', 'Total_Dissolved_Phosphorus', 'Ammonia', 'Nitrate_Nitrite']):
        cc = np.corrcoef([nutrient_monthly[:,i+1].astype(float), runoff_shortened['flow']])
        coef = cc[0,1]
        rho, p = kendalltau(nutrient_monthly[:,i+1].astype(float), runoff_shortened['flow'])
        
        # save in matrix!
        CC_all[i,j] = coef
        Kendall_all[i,j] = rho
        
print(CC_all)
print(Kendall_all)

[[ 0.55171746  0.05101281  0.42698623 -0.06435019  0.71538914  0.42565129
   0.03099001]
 [ 0.74625689  0.38591678  0.80995668  0.67482507  0.77028754  0.70080655
   0.16274738]
 [ 0.43621572  0.33513288  0.61280153  0.53963714  0.34177486  0.35913761
   0.34263741]
 [ 0.84486963  0.41300529  0.76362528  0.74787927  0.67312942  0.83938023
   0.73951505]]
[[ 0.25581395 -0.06341463  0.22702104  0.02325581  0.2         0.10107527
  -0.08965517]
 [ 0.35326689  0.04878049  0.45736434  0.33554817  0.57894737  0.2516129
   0.00689655]
 [ 0.24473976  0.10243902  0.38648948  0.35105205  0.41052632 -0.01075269
   0.00689655]
 [ 0.52159468  0.35265412  0.5503876   0.48837209  0.46315789  0.54251915
   0.48505747]]


# make the color plot!
fig, (ax1, ax2) = plt.subplots(2, 1, figsize = (4.5, 6))
plt.subplots_adjust(hspace=0.6)

# Create the "pcolor" object while plotting it to the axes.
# Use our pcolor object to create a colorbar.
pcol = ax1.pcolor(CC_all, cmap='PRGn_r', vmin=-1., vmax=1.)
cbar = plt.colorbar(pcol, ax=ax1)
cbar.set_label('Coefficient')

# Set some better ticks and tick labels!
xlabels = ['WE2', 'WE4', 'WE6', 'WE8', 'WE9', 'WE12', 'WE13']
ylabels = ['TP', 'TDP', 'A', 'N/N']

# Now, set our ticks and labels. Everything lines up!
ax1.set_xticks(np.arange(0, 7, 1) + .5, xlabels)
ax1.set_yticks(np.arange(0, 4, 1) + .5, ylabels)

# more formatting
ax1.tick_params(top=True, labeltop=True, bottom=False, labelbottom=False, pad=5)
ax1.set_title('Correlation Coefficient: Runoff vs. All Site Nutrients!', pad=10)

# Do the same for the Spearman rank coef. matrix:
pcol = ax2.pcolor(Kendall_all, cmap='PRGn_r', vmin=-1., vmax=1.)
cbar = plt.colorbar(pcol, ax=ax2)
cbar.set_label(r"Kendall's $\tau$")

# Set some better ticks and tick labels!
xlabels = ['WE2', 'WE4', 'WE6', 'WE8', 'WE9', 'WE12', 'WE13']
ylabels = ['TP', 'TDP', 'A', 'N/N']

# Now, set our ticks and labels. Everything lines up!
ax2.set_xticks(np.arange(0, 7, 1) + .5, xlabels)
ax2.set_yticks(np.arange(0, 4, 1) + .5, ylabels)

# more formatting
ax2.tick_params(top=True, labeltop=True, bottom=False, labelbottom=False, pad=5)
ax2.set_title("Kendall's \N{GREEK SMALL LETTER TAU}: Runoff vs. All Site Nutrients!", pad=10)

plt.show()

Final Project: Lake Erie & Harmful Algal Blooms¶

Sara Chojnowski, Merjem Memic, Alicia Schmidt, Jessica Zhang¶

EARTH/SPACE/CLIMATE 323¶

April 19, 2024¶

Purpose and Introduction¶

Primary Research Questions¶

1. What are the primary trends and variations in Lake Erie nutrient levels between 2012 to 2018? Observe and characterize the behavior in nutrient concentrations over time, taking note of seasonal variations, relationships between variables, and HABs events.¶

2. What is the correlation between upstream nutrient release (agricultural runoff & wastewater discharge) and downstream nutrient concentrations? Further, how does this coincide with recorded HABS events?¶

Data¶

Approach and Results¶

PART 1: Characterizing Nutrient Data & HABS Events¶

Question 1¶

Question 2¶

Question 3¶

Question 4¶

Question 5¶

PART 2: Correlating Upstream and Downstream Nutrients¶

Question 6¶

Question 7¶

Question 8¶

Conclusions¶

References¶